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Abstract. We study the Glauber dynamics of Ising spin models with random 
bonds, on finitely connected random graphs. We generalize a recent dynamical 
replica theory with which to predict the evolution of the joint spin-field 
distribution, to include random graphs with arbitrary degree distributions. The 
theory is applied to Ising ferromagnets on randomly diluted Bethe lattices, where 
we study the evolution of the magnetization and the internal energy. It predicts 
a prominent slowing down of the flow in the Griffiths phase, it suggests a further 
dynamical transition at lower temperatures within the Griffiths phase, and it is 
verified quantitatively by the results of Monte Carlo simulations. 



1. Introduction 

Finitely connected (FC) spin systems were introduced more than 20 years ago by 
Viana and Bray pQ as more realistic alternatives to infinite range (IR) models of spin- 
glasses [H |3] . In finitely connected systems the spins are placed on the vertices of a 
random graph, and interact only when their vertices are connected; the number of 
connections per spin remains finite (on average), even in the thermodynamic limit. 
This definition endows finitely connected spin models with a geometry (e.g. vertex 
neighborhood), a crucial feature also of finite dimensional (FD) spin systems, that 
was absent from infinite range models. Yet, in contrast to FD spin systems which 
arc notoriously difficult to solve, FC models are still of a mean field nature and 
can therefore be studied analytically using methods from the statistical mechanics 
of disordered systems. This property reflects the absence of short loops: in finitely 
connected spin systems loop lengths are typically of order log(iV), so that the spins live 
in environments which are locally tree-like, unlike spins in finite dimensional systems, 
and short-range frustration cannot occur. As a result of their analytical accessability 
the equilibrium properties of finitely connected spin systems are now understood quite 
well [H O [6l 13 [H [9] . The mathematical and numerical techniques which originated 
from these equilibrium papers were, in turn, generalized and applied in subsequent 
dynamical studies [IHl [IT1 [121 [13l [III [15l [16l [13 [18] . 

One of the properties shared by finitely connected and finite dimensional spin 
systems is the presence of Griffiths singularities [19 . In his seminal paper Griffiths 
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showed that in the diluted Ising ferromagnet, where either sites or bonds of a classical 
lattice are removed with some probability 1— p, the magnetization is a non-analytical 
function of the external field for a range of temperatures T c (p) < T < T c (l), 
where T c {p) and T c (l) are the critical temperatures marking the P-^F transition 
of the diluted and undiluted systems, respectively. The system is in a conventional 
paramagnetic state only for temperatures above T c (l), where T c (l) could be infinite 
[20] . The temperature interval T c (p) < T < T c (l) over which these singularities occur 
is called the Griffiths phase [21]. This peculiar behavior of the magnetization [22] 
and other thermodynamic functions is understood to be caused by the presence in the 
randomly diluted system of large undiluted spatial regions (or clusters) of the lattice. 
In the Griffiths phase these clusters are in an ordered magnetic state, although the 
system is globally paramagnetic. The Griffiths singularities are not always strong 
[22] [24] [25] [20] and often difficult to observe experimentally [24] , nevertheless this is 
possible with modern sampling techniques [21)] . 

In contrast to statics, the effects of large undiluted clusters on the dynamic 
properties of diluted spin systems are more drastic. The dynamics in such clusters 
is very slow because it requires reversing spins coherently in the entire cluster. In 
FD spin systems this results in non-exponential decay of the spin autocorrelation 
and magnetization functions in the entire Griffiths region [21] [27] [28] [29] [30] 131] . 
The latter studies concentrated mainly on the derivation of bounds for the spin 
autocorrelation function, at large times, with subsequent verification by Monte Carlo 
(MC) simulations. The dynamic properties of the Griffiths phase in FC spin systems 
remain, to the best of our knowledge and that of others [32], largely unexplored. 

In this paper we generalize recent results obtained within the framework of 
dynamic replica theory (DRT) 17J to include random graphs with arbitrary vertex 
degree distributions, and apply the generalized theory to the dynamics of diluted 
ferromagnets in the Griffiths phase. Our paper is organized as follows. In section[2]we 
define our finitely connected spin model and its dynamical equations. In section [3] we 
close the macroscopic dynamical laws using the standard assumptions and procedures 
of DRT. From these closed laws we recover known results of equilibrium statistical 
mechanics as stationary solutions in section 13.21 as a test. The replica-symmetry 
assumption allows us to take the replica limit n — > in section 13.31 The resulting 
dynamical theory is applied to the Glauber dynamics of diluted Ising ferromagnet in 
section |H We close with a summary and discussion of our results. 

2. Model definitions and dynamic equations 

We consider a system of N Ising spins ct% <E { — 1,1}, which are placed on the vertices 
of a finitely connected random graph. Spins interact only when they are connected. 
Their microscopic dynamics is described by a master equation for the evolution of the 
microscopic state probability in continuous time: 



in which <x = (ci, . . . , o"at), Fi denotes the spin-flip operator defined via FiVL(cr) = 
r2(eri, . . . , — crj, . . . , ctat), and the quantities Wi{cr) are the Glauber transition rates 




(1) 



w i( <T ) = - Cjtanh[/3/ii(er)]] 



(2) 



\ See 1231 for a model example where very strong Griffiths effects arc found. 
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with the local fields 

hi(<r) = c ijJij°~j + Q- (3) 

The parameters (3 = T^ 1 and 9 define the inverse temperature and a uniform 
external field, respectively. The random interactions {cijJij} are symmetric, viz. 
CijJij = CjiJji, and are regarded as a quenched disorder. The interaction strengths 
Jij are independent random variables, drawn from a probability distribution P(J). 
The random variables Cy S {0, 1} are the entries of an adjacency matrix, with zeroes 
on the main diagonal, defining the random interaction graph. The symmetry of the 
interactions ensures that the process ([1]) evolves towards equilibrium, characterized 
by the Boltzmann measure Poo(f) ~ exp[—/3H(tr)], with the Hamiltonian 

H(tr) = -^(TiCijJijO-j - 9 J~] <Jj . (4) 

i<j i 

In this paper we consider FC random graphs where the vertex degrees {ki}, with 
ki = J^j^i c iji are drawn randomly and independently from an arbitrary probability 
distribution P c (k) over the non-negative integers, with finite first moment c = 
^2uP c (k)k. The probability of finding an adjacency matrix c = {c^} in this random 
graph ensemble, constrained by the vertex degrees {ki}, is given by 

P(c\{ki}) ^ 11/ « , s II > >: , < (5) 

i<j i 

where Z is a normalization constant, and 

The presence of p c (c«) in the definition (O is mathematically convenient in solving 
the model, but not essential; it can be transformed away in leading order in N. 

We avoid the impossible task of solving the 2 N equations fT} directly, and consider 
an alternative description of the dynamics in terms of macroscopic observables. In 
particular, we consider the evolution in time of the joint spin- field distribution [33], 
which is given by 

D(s,h;a) = ^J2 6 ^ S l h - h ^- ( ? ) 

i 

In finitely connected models equipped with the dynamics ([1]), the macroscopic 
distribution J7J will evolve deterministically for N — > oo, according to a macroscopic 
dynamical equation [17] of the form 

5 1 1 

— D(s, h) = - [1 + s tanh[/?/i]] D(-s, h)- -[l-s tanh[/3/i]] D(s, h) 

1 X -, 



(8) 



Mi<j: p c (cij) = — 5^,1 + (1 - ^)<5 c «,o • (6) 



2 C 2^ / dh '^ - s '^MP h ']\ A \^ s '\ h : h '\ s '] 

s' 

J ^'[l-s'taxih[f3h']}A[s,s';h,h';0] 

s' 

with a spin variable s G { — 1,1}, and a field h G K. The dynamical equation ([8]) is 
written in terms of time- dependent kernels D and A, which are defined as follows 

D(s,h) = ^£<W[/»-M*)]>A-t (9) 
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(10) 



e.t' 



A U '[s,s';h, h';§] = (^ s > <ai S s<ae/ 5[h'-h i (a)]5[h-he(cr)+2J u >s]j ^ 

with s' 6 { — 1,1}, h! <E R and s € {0, s'}. In these expressions, the sub-shell average 

_ E<r PtWto FU « h) - D(s, h; a)} 



</(*)>; 



(11) 



/D;t E < T'Pt(»')n.fc*^(*.^)--D(*.^^)] 

is written in terms of the macroscopic distribution ([7]) which acts as a constraint 
on micro-state^], and the microscopic probability distribution p t {er). The kernel A 
is positive semi-definite, and normalized for N — > oo; it defines the joint spin-field 
probability distribution of connected sites. Equation is exact for large N, but not 
yet closed due to the presence of microscopic probability pt(er) in (jTTJ) . 

3. Dynamical replica analysis 

3.1. Closure and disorder averaging 

In order to solve equation (jHJl we have to compute the kernel ([TO)) . This latter kernel 
is dependent on the disorder {c.y J^} and the microscopic state probability pt(cr). To 
compute A we make the usual assumptions of the dynamic replica method |331 117) : 
(i) the observables {D(s, h; <x)} are assumed to be self-averaging with respect to the 
disorder at any time, and (ii) the microscopic probability pt{cr) is taken to depend on 
er only through {D(s, h; tr)}. The self-averaging assumption leads us to 



A[s,s';h 7 h';s] = (^y^cu'Au'[s,s';h,h';s\) 



(12) 



The subsequent elimination of the microscopic probability Pt{&) from the 
above, followed by the elimination of the fraction via the replica identity 
£<r $(tr)W(*)/ Y, a , W(ct') = lim„^ . . . <Z>(^) n£ =1 W(<r a ) allow us to 
perform the disorder averages in the term cu'Aw of equation (|10[) (see |Appendix A| 
for details), yielding 

{cu>Au> [s, s'; /i, ft,'; s]){ 0w j y } 



E--E /n[ d ^ d ^ ^kHf\ 

x H 6 [d(t, h)-±J2 6 ^? 5 i h - H ° 

rha i 

x 5 s ,^8 s ^8[ti - H}] (c w 5[h - Hj, + 2J, 



-iE„ Khi(<T a ) 



{C-i j Jij } 



1 c 
ZN 



E-E/n 



dfffdft? 
2^ 



§ Here, to simplify notation, we skip explicit mentioning of the intermediate discretization of the 
fields h in J7), which is formally required 1 171 ■ 
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n 



dk 



2tt 



-i{ke+k e /} 



x exp 



7^E{ f dJP ( J ) c- iJ ^- {' h f^+ h 7<}- i ^+^ -I'j + 0(N°)Y (13) 



In the derivation of the above result we have used the integral representation of unity 
l=fl[dH?8[H?-h i {<r a )] (14) 



and the integral representation of the Kronecker (5-functions 



(15) 



In order to disentangle the N degrees of freedom in equation (flU)) we next define a 
replica density function 

P(<7, h, h {<7 4 }, {hi}, {h}) = ^ £ - ~ *i] ( 16 ) 

i 

where <r = (ai, . . . er„), <Tj = [a\, . . . ex") (similarly for the replicated vectors h, etc.), 
via insertion into equation (|13[) of the following (5-functional unity representation 

1=1 J] dP(o-,A,fe)*[p(o-,A,ft)-P(o-,h,^{o-i},{Ai},{fe})] (17) 

which gives, with the short-hands (g(J))j = fd J P(J)g(J) and x • y = J2 a x a y a , 
(cm Aw [8, a'; h, h'; s]) {CijJtj} (18) 



1 c 

Z7V 



n 



dO a (r, /Q 
2-k/N 



n 

cr,h.k 



dP(<T,h,k)dP(tr,h,k) 

27r/iV 



exp [iN^2 D a (T,h)D(T,h) +iN^2 J dhdkP{(r,h,k)P(cr,h,k) 

T,h,a <7 

-cN [ dhdhdkdk'P{er, h, k)P(a', h ', k') 

2 cr criJ 



K ^ c -iJ[h-cr'+h ,(T]-i[k+k'] - 1^ + 0(N°) 



ee /n 



rd-ffjd/i, 
I 2tt 



77 



e^]ex P [i£M^-0} 



exp [-i ^ DaCr.^^^-tfffc-flfl-i^PCo-i.hj.fcj) 



Inserting the above result into the sum (jlOp . followed by further manipulations (see 
| Appendix B| for details), leads us to the path integral 

N 

<{dPdPdD} e »mP,P,i>}]+o(N°) 



1 



A\s, s : h, h: s] = lim lim — 



1 

27 
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r h'~>n rr n-i J 



k,k'>0 

xM[H, h,cr\k-l, 9} M[H', h , a'\k'-l, ( 

xS^SisJlh'-H!] (s[h-H[+2Js] e^Jlh.tr'+h'.aA (ig) 



[dHdH'dhdh M[H,h,a\k, 6} M[H',h', a'\k',6] '+ 0(JV 4 ) 
a. a' •* 

where we use the following definitions: 
*[{P,P,£>}] = i ^ AhD a (T,h)D{r,h) +ij^ J dhdkP(o-,h,k)P(o-,h,k) 

r,h,a (T 

+ -cJ2 [dhdhdkdk'P(a,h^k)P(a\h\k')(e- u ^ 
cr.cr' 

+ ^P c (fc)log^ j dHdhM{H,h,a-\k,e] (20) 
fe>0 cr J 

and 

dfc e -ifcm M[H, h, a\k, k, 6} (21) 

M\H h erlfc k 9} = — e ^.{H -6}-\Y. T „ a AhD a (r,h)s^ C7a s[h-H a ]+ikk-iP((T,h,k) 
' ' 2tt 

with m S Z. Finally, we change the order of the limits N — > oo and n — » in (UHJ) 
and, with the help of the normalization identity ^ s s , J" dM/i'A[s, s'; /i, ft,'; s] — 1, we 
compute (HHJ) by steepest descent, which gives 

A[s, s; h, h'; s] = lim — ^ P c (k)P c (k') ^ [ dHdH'dhdh 
A k,k' a, a 1 •* 

x M[H,h,tr\h-l,6] M[H', h, cr'\k'-l, 9} 

x S sl , ai S s ^5[h' - Hi](6[h -H[ + 2Js]e- iJ ^- <T '+ h '- (T ^ j 

x[Yl [dHdH'dhdh M[H,h,a\k, 6] M[H',ti ,tr'\k',6] 1 (22) 

where Za is a constant that ensures the proper normalization of A. The order 
parameters {P,P,D} are determined by extremization of the functional 4" in (|20|) . 
which leads us to four functional saddle-point equations 

^V^^^/^M^,^,^^^-^) 

P( ff! M) = EfiW j^glA^lM (24 ) 
P(er,h,k) = icQ(a;h,k) (25) 
Q(er, M)=£ Jdh'dk'P{a', h, k')^ e -iAh.<r'+h' .a-]-i[k+k'] _ ^ (26) 
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The relations (|25I26|) allow us to relate the order parameter P(cr, h, fc) to its 
conjugate P(cr,h,k), and thereby remove the latter from the function M (|2~lj) . 
Furthermore, assuming that the function D a (s,h) is well behaved in the sense that 
J\ Ah D a (s, h)g(h) — > Jdh D a (s, h)g(h) for Ah — » 0, we have 

M[H, h, a\k, k,0] = — ffi-{**-0}-iX a D«(*M+ikk+cQ«T,h;k) (27) 
in definition (|2Tj) . 

The conjugate parameters D a (a,h) and P(cr,h,k) in our replica theory play 
the role of Lagrange multipliers enforcing the normalization of the joint spin-field 
distribution D(a,h) and of the replica density function P(<r,h, fc). The physical 
meaning of the density P(<x, h, fc) is not yet clear due to the presence of the vector h 
and the parameter fc. However, we note that in our theory only the Fourier transforms 

Jdh e~ ix h J\ dfc e- ik P(cr, h, fc) of this function are relevant, where x £ K™. 
3.2. Equilibrium 

In this section we show that the equilibrium solution of the model (j4|) is also a 
stationary solution of our dynamic equation (jSJ. This is done in two steps. First, 
we show that the equilibrium replica theory of the model under study is a special case 
of our dynamical replica theory. In order to do this, we make an ansatz as in [33] : 

and evaluate the Fourier transform of the replica density 



dhe~ ix fl dfc e- ikm P((T, ft, fc) (29) 

J —7T 

for x £ R™ and m G Z. Using the saddle-point equation (|24|) for the order parameter 
function P{cr 1 h, fc), combined with the Fourier transform of the function M (see 
| Appendix C| for details), 

dfc e~ ii:m M[H 7 h 7 cr\k,k,8} 

fc-r 



g— Cgfc— m 



n 



., iJL J2 I dh e dJeP(J e ) [ dk e P(<T i ,h i ,k i )e- i ' ke e- iJ t fl '- (T 

\ l = l \_ (Te J J -TV 

(27r) n s\H-J2 Jt<rt-0-x] e- 1 ^^-^-.^-) (30) 



leads us to the desired result for 

dhe~ ix fl [ dfc e- iiem P{<r, h, fc) (31) 



fc! 



k>7 
k-i 

x 



M fc (fc-m)! 



n 



I 'dhidJiP(Jt) f dk e P( ( r e ,h e ,k e )e- iki e- iJifl '- (T 
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Solving equation (|3ip for m = 1 yields a very useful equality, 



E / dhidJiP(J t ) / dk e P((Ti,h e ,k e )e- [ke e- [ 



(32) 



-iX.fl 



dk e _i *P(«r, £, £) = / dfc e -iJ P(o-, fe) e^ CT :c 



d/i e 



which allows us to compute the integrals over {hi} in pip , giving 



dft. e 



ix.h 



dk e- ikm P{a-,h,k) 



= E 



A:! 



X: (T ,P(cr')/dJP(J) c f Jcr cr ' 



k—r 



(33) 



(34) 



e /3<r.0 e i/3cr.a; 



(k — ra)\ 



E<r< E<T>»P(<7'")JdJP(J)e^"v 



^/3<T".6 



with the shorthand P(tr) = J* dk e~ lfc P(er, k). Now for x = (0, ... , 0) and m = 1 
our equation (|34[) acquires the form 

fc-i 



w = E 



Pc(fc)fc 



'E tr ,P(<r')fdJP(J) e P J(T(T ' 



k>l 



Der" Y,*»>H<r" , )J*JP(J)eP J ""-< r '" 



P /3(T".0 



(35) 



which is exactly the order parameter equation as found in equilibrium [34] . 

The second part of our proof consists in showing that the ansatz (|28|) also leads 
to a stationary solution of our present dynamic equation (|8|). For this purpose we 
compute the saddle-point equations for the joint spin-field probability distributions 
l]22p and (|2"5)h given our ansatz. The result of this calculation (see |Appcndix D] for 
details) allows us to write these two equations in the form 



D(s,h) 
A[s, s ; h, h ; s] 



- o 0sh 



$[h] 



fls(h+2JS)+f3s'h' -f3Jss' 



A[h+2Js-Js';h'-Js] 



(36) 
(37) 



respectively, where all complicated terms dependent on the replicas are contained in 
the two functions 4> and A (which are defined in | Appendix D[ ). Inserting (j3l)|) and (f3"7| 
into the right-hand side of the dynamic equation |8[). followed by further manipulations 
(see | Appendix E ), leads us to the equality Jj£)(s, h) = for all s G { — 1, 1} and Ii6l, 
as claimed. Thus, the equilibrium solution of the model ^ indeed defines a stationary 
state of the dynamics ijHJl . 



3.3. Replica symmetry 

In order to take the n — ► limit in our equations (|22p ~ (|24p we assume replica symmetry 
(RS). For the conjugate order parameters P Q (s, H), which are depend only on a single 
replica index and expected to be imaginary, this translates into 

D a (s,H) = ilogd(s,H). (38) 

The replica density P(c, h, k) depends on one discrete vector er and one continuous 
vector h in replica space. The parameter k is a scalar variable coupled to the vertex 
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degree k 7 which is a random variable. Replica symmetry demands that the order 
parameter P(<r, h, k) is invariant under permutation of the replica indices, for any 
value of k, which implies [H [35] that it is of the general form 



/n 
{dP} W[{P};k] J] P{a a ,h a ) 
a=l 



(39) 



where ^[{Pjjfc] is a functional distribution, which must be normalized according 
to J{dP} dk W[{P}; k] = 1. It turns out that also the Fourier transform 

die W[{P};k] e~ lk of this functional distribution is normalizeojjj], which is very 
convenient for our further calculations. 

The RS ansatz (138139)) allows us to take the replica limit n — ► in equations 
P2 )l -I(2g p . The Fourier transform J^dk e~ ikm M[H, h, a\k, k, 9], where m 6 Z, is the 
main ingredient of these equations. We can easily compute its RS version using result 
(|3"U)) . leading to 



dk e- ikm M R s\H 7 h 7 cr\k,k,, 



(40) 



£=1 

A:— m 



(& — m)! 



Now we can use (|4"0")) and the saddle-point equation (|2"4")) to solve for the functional 
distribution W[{P};fc]. However, it is clear from (|4"0"1) that all the equations of our 
theory are dependent only on dk W[{P}; k] e~ lkm , rather than on the distribution 
^[{P}; k] itself. Thus for m G Z we define 



W[{P}\m] 



dk W[{P};k] e 



-ikm 



(41) 



and compute this object (see | Appendix F| , which leads us to the equation 

„k—m 



W[{P}\m] = Pc(k) 



k>m 



P(a, h)- 



(k — m)\ 



J] [dJ t P(Jt){dP e }W[{P t }\l]] 



(42) 



Z[{P!,...,P fe _ m }] 

where m G {0, 1}, and Z[. . .] is a normalization constant given by 



k — m 



Z[{Pi, ■ 



1 Pk — 711 



53 / MiPtfaht) e~ ijJl ^' 



d(a',J2Jio- e + 9) 

i 

(43) 



| This can be shown by substituting J39H into the function \t (120 ft , followed by expanding this 
function for small n. The desired result f{dP} dk W[{P}; k] e~ lk = 1 then follows from solving 
the saddle-point equations for the O(n ) part of 
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It is easy to verify that for m € {0, 1} the functional distribution (|4"!?|) is normalized 
for any vertex degree distribution P c (k), provided the latter satisfies J2 k P c (k)k = c. 

Next we compute the kernels (|22p and (|2"5|) under the RS ansatz. This is done by 
using (|4"0"1) in both, followed by the replica limit (see |Appendix F[ ) , giving 

D(a,h) = J2 p d k ) I II [*JtP{Jt) {dPt} W[{P t }\l]] (44) 



fc>0 











d(a,J2e Jeo-£+6) 



and 

A[s, s'; h, h'; s) 



E 

k,k'>l 
k-1 



P c (k)k P c (k')k' 



J dJ P(J) 



(45) 



[ Y[[dJ e P(Jt){dP e }W[{Pe}\l]}] [ / II [dJ' r P{J' r ){dQ r }W[{Q r }\l]] 

•* 1=1 J r=l 



k'-l 



k-1 



En 



a,a' 1=1 
k'-l 



/ dhtPt(ai,he) e" 



r=l 



^ / d/t r Q r (<r r ,/i r ) e - iJ >"' T ' 



cZ(s, /i + 2Js) 



S ,I7 U S,(J 



5[ti- -hot -0- Ja']S[h - ^ 4°> - - Ja + 2JI 



k-1 . 

Ell [E dhePda^hi) e- iJ « h * a ]d(a,J2Ji<re + + J°') 



K a,a' 1=1 o t 
k'-l 



II [ E / dh r Qr(<J r , K) e - iJ ^" CT '] d{a', ^2 J' r cj r + 9 + Ja) \ 

r=l oy r ) 

where s £ {0, s'}. Equations (|42 |) - (|45|) constitute the final analytic results of the 
RS theory in this section. The results of a similar dynamical study [17j . which was 
carried out for Poissonian graphs only, are easily recovered from the present more 
general equations, by using the equality 

E P <W c ~ m = E P ^ k ) «*- ( 46 ) 

k>m ^ m '' k>0 



(which holds for all m € {0, 1, . . . , } for the Poisson vertex degree distribution, i.e. 
when P c (k) = e~ c c k /kl), throughout formulae (|47]) - (f4"S]) . 

The solution of our dynamic equation ([8]) requires the computation of the kernel 
(|45p at every instance of time t. In order to compute this kernel we have to solve the 
saddle-point equations (JUJ) and (|4"4")) for the functional distribution W and the function 
c?(s,/i), given the instantaneous values of the joint spin-field distribution D{s 1 h) at 
time t. These equations are integro-functional, and analytic solution is generally ruled 
out. However, we can solve them numerically [17) using the population dynamics [7] 
algorithm. In order to apply this latter numerical method efficiently we may transform 
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W — > W and P(o~\x ) -> jdhP{a,h)e- ihx in equations (|g2 ]) -([15 )) . according 
W / [{f'}|l]= / {d^}VK[{P}|l]J^(5 P{a\x)- J dhP{a, h) e 4 * 



(47) 



Upon substitution of (|42p into (j4"T|) we can easily derive the functional relation for 
(|47p . which is given by 

. fc-i 



W^[{P}|1] = E / II {d^W){dA}W[{A}|l]} (48) 

WjZI {E<t £ A(^l-Tieo-)} ^(^Efci 1 Jio-e+O+x) 



k>l 



X 



P(a\x) 



Z[{P u ...,P k ^}} 



where ■ ■ ■ , Ph-i}] = £ CT Ilt^E^ A(^|J^)}rf(a, E^i ^ + 0). The 

normalization of W is seen to be built into this relation, however the functional 
arguments P(o~\x) are only normalized for x = 0. 

4. Dynamics of the diluted Ising ferromagnet in the Griffiths phase 

As an explicit application of the theory derived in previous sections, we now study 
the Glauber dynamics of the diluted Ising ferromagnet on the Bethe lattice. 

4-1- The model and its equilibrium properties 

We consider a model of an Ising ferromagnet characterized by the following 
Hamiltonian: 

(ij) * 

The sum is over all links of the Bethe lattice with connectivity k. The bonds Jij are 
random and statistically independent: Jjj = J with probability p and = with 
probability 1 — p. The lattice contains only finite size clusters for p < p c , where p c is 
the percolation threshold given by p c = l/(k — 1) for the Bethe lattice [36], whereas 
the giant cluster appears for p > p c . The density of the finite clusters of bond-size n 
is also known for the present model [36j . and asymptotically given by 

W n (p, k) ~ n ~ie~ nA ^ (n -> oo) (50) 

where 



A(p, k) = In 



(fc-2) 



fc-2 



(51) 



(k- l)^ 1 p(l -p) fe " 2 . 

For p = p c we have A(p, fc) = and the asymptotic form (f50|) of the density W n (p, k) is 
independent of fc. The model (|4^| has paramagnetic and ferromagnetic phases, which 
are separated by the critical boundary [3"T] 

T c (p) = J/tmh-\p c /p). (52) 

The critical temperature of the undiluted Ising ferromagnet on the Bethe lattice is 
simply T c (l). Thus the Griffith phase of the model (|49|) is given by the range of 

% Here we assume that for x S K the Fourier transforms JdhP(a, h)e~ lhx are real-valued, which is 
certainly true in equilibrium (see section [3.21 1 and is a self-consistent assumption for any time. 
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temperatures T c (p) < T < T c (l). The magnetization in the Griffiths phase and in the 
paramagnetic phase (i.e. for T > T c (l)) vanishes, and without an external field (i.e. 
for 9 = 0) the internal energy is given by 

(ff(tr)) = -ipfctanh(J/T) (53) 

where the angular brackets define a thermal average (expression ([53]) is easily derived 
from the free energy in [33] ). The presence of Griffiths singularities in the low 
temperature part of Griffiths region was demonstrated in [37j by studying the density 
of Yang-Lee zeroes [351 13S]- Moreover, the authors of [37] obtained an exact expansion 
for the cluster magnetizations, which was used in arguments by Harris [22j for the site- 
diluted version of this problem, within the cavity approach. The presence of rare large 
clusters in the diluted Bethe lattice is responsible for the Griffiths effects in this model 
[37] . This singularity, however, is very weak (<~ e - const /l e l) an d j s difficult to observe in 
equilibrium [24] . In this paper we consider the Glauber dynamics of the diluted Ising 
ferromagnet (|49[) in the paramagnetic and Griffiths phases. To connect our dynamical 
theory, which was developed for random graphs parameterized by an arbitrary vertex 
degree distribution, with the equilibrium studies of this problem as carried out for 
Bethe lattices, we note that in the infinite system size limit N — > oo the random 
regular graphs defined by P c (k) = 5}- tC asymptotically approach Bethe lattices [3D]. 



4-2. Equations of the DRT for random regular graphs with dilution 

We can derive the order parameter equations (|42I44I45P for the diluted Ising 
ferromagnet simply by inserting into these three general equations the special choices 
P c (k) = 5k, c and P(Je) = P 5(Jg — J) + q 5(Jg), where p e [0, 1] and q = I— p. Equation 
l|42[) for the order parameter function W is then simplified by summation over k. If 
we also replace c — ► k (since the non-diluted graph is regular), this leads us to 

k— 1 

W[{P}\1}= Pin,..-,Tk-i) fl[[{dPt} W[{P t }\l]] (54) 



i=i 



a J) 



P(a,h) 



JdH d(a,H)e 



ih\H- 



llLi 1 EaJdhePeiae/he)^ 



-iJrg \hi<T-\-h<Tg\ 



where we define the probability function 

P(n,...,T k ) =pELi^.igfe-E?=i^.i (55) 

with the binary variable r G {1,0}. We note that in equation (|54]1 the terms with 
ti = do not contribute, since the Pi(ai^ he) are normalized by definition. Finally, we 
transform W[{P}|1] — > W[{P}|1] in equation ((54]) . according to the definition 



w[{p}\i}= / {dPi^piiiin^ 



P{a\a') 



dhP(a,h)e 



-ihja' 



(56) 



where a, a' G { — 1, 1}, which leads us to an equation for the functional distribution of 
Fourier transforms: 



W 



k — 1 

{P}|1]= ]T P(r 1 ,...,7 )b _ 1 )/n[{dP / } W[{P t }\l] 



(57) 



ri,...,r fc _ 1 
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We can in fact get rid of the Tg variables entirely, which gives us an alternative 
representation of the equation above 



W[{P}\1] = £ Bk-x&) [ fl [{ d ^} W[{P t }\l] 

fe'=0 t=i 



(58) 



PW) 



Till {J2 ae PtfrlW)} d(a, J Eli <rt+6+ J°') 



Eer" ^=1 {E„ AfoK)} d(a", JEtl ' 

where Bk-i(k') is the binomial distribution 



fc-1 
fc' 



p k' q k-i-k'_ 



(59) 



This result reflects the fact that the distribution of the vertex degrees in the random 
regular graph of degree k with the bond-dilution is indeed the binomial Bk(k'). The 
fields ((3]) for the model (|49|) take the values Jrt + 8, where n S {— A;, . . . , fc}, allowing 
us to write the joint spin-field probability distributions and in the form 



D(s,h)= P(s,n) S(h- Jn-9) 



(60) 



fe-i fe-i 

A[s, s'; ft, A'; s] = ^ ^ /a[s, s'; n, n'\r] 5[ti - Jri - - Jts] 

n— — k-\-l n' — — fe+1 

x 5 [ft + 2 Jr5 -Jn-9- Jts'] 
where (. . .) T = J2 T P( T )i with P(t) defined in (f55|) . and 

P(s,n) = X! P(n,...,T k ) JU^dPi) W[{P t }\l] 



(61) 



(62) 





^|r £ s) 




E CT n^=i 




d(cr,JX^T^ + 0) 



A[«,a / ;n,n'|T]= ]T P(n, . . . , t*_i) /jj [{dft}W[{P/}|l] 

n,— ,r«i-i £=i 

x ^ P(r{,...,rU) /if [{d&W{0r}|l] 

fe-i r 

X II XI Mains') d(s', Jri + + Jts) S n ,.^-i 
fe-i r 

X IT ^ZQr{o- r \T' r s) d(s,Jn + 8 + Jts') S n .^ k -i T 



(63) 
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fe-i 



En 



a.a' 1=1 

fc-1 

X 
r=l 



n 



. Of, 

YjQrio-rVW) 



fc-1 

d(i7, j y2 T z a t + ° + Jt<jI ) 
fe-i 

d(cr' , J^~^ r^cr,. + + Jrer) 



where in deriving probability distributions over the integer fields (|62p and (|63|) we 
followed the steps leading to (p)Tj) . It is easy to show, using equation {57]), that the 
distribution P(s, n) is the marginal of ^ T A[s, s'; n, ti'\t]P(t). The simplified form of 
the probability distributions ([HO)) and (|6"Tj) allows us to reduce our dynamic equation 
(8j to a system of ordinary differential equations (see |Appcndix G| for details) 

d P t {s,n) = \ [1+s tanh[/3 Jn + (36}} P t (-s,n) - ^ [1 - s tanh[/3 Jn + [36]} P t (s, ri) 



dt 



2 L ' JJ v ' " 2 

fc-i 

+ A t [s, l;n + l,n'|l] -[1 - tanh[/3J(n'+ s) 



n'=-fc+l 
fe-1 



+ ^ A t [s, -l;n - l,n'|l] ^[1 + tanh[/3J(n'+ s) - 

n' = -fe+l 

fe-1 _^ 

-pk Y A t [s, -l;n+l,n'\l} -[1 +tanh[^J(n'+ s) • 

n' = -fe+l 

fe-1 . 

-pfc ^ A t [s,l;n-l,n'\l] -[1 -t&nh[f3J(n'+ s) + [38]}. (64) 

n' = -fe+l 

Here n s {—A:, . . . , /c}, and A t [s, s'; n, n'|l] = for n, n' ^ {— k + 1, . . . , k — 1}, leading 
to four boundary equations. The equations of the dynamical replica theory (|5TI62I63|) 
and (|64j) are now cast into a form which allows us to solve them numerically. 



4-3. Numerical results 

Here we use the analytic results of the previous section to study the dynamics in 
the Griffiths phase of the diluted Ising ferromagnet (|49|) . We solve the dynamical 
equation (|64[) for the probability distribution P t (s,n) numerically, given the initial 
values (see |Appendix H[ ) and given the boundary conditions of this equation, using 
Euler's forward iteration method. At each iteration step of this method we solve 
equations (|57p and (|^|) , using a population dynamics algorithm (see |Appendix I| ) , 
for the distribution W and the function d. The result is then used to compute the 
probability distribution (|63[) and to iterate the discrete version of the dynamic equation 
(|64| over the next time step t — > t + At. In order to assess the quality of our dynamic 
theory we compare results of our numerical solutions of (|64p with the results of Monte 
Carlo (MC) simulations. In each simulation we generate a random regular graph of 
degree k with N vertices using the algorithm of Steger and Wormald |41j . We then 
remove each of the edges from this graph with probability 1 — p, so that on average 
only pkN/2 edges remain in the resulting diluted random graph. Finally, we perform 
MC simulations of the ferromagnetic Ising model defined on the diluted random graph 
P5j) using conventional Glauber dynamics. 

The evolution in time of the magnetization and the energy per spin, as obtained 
firstly in the numeric solution of theory and secondly in the MC simulations, is 
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P(-l,n) 



Figure 1. Left: evolution of the magnetization and energy per spin for k = 3, 
p = 0.2, J = 1 and = 0. The temperatures are T = 3 (bottom lines) and 
T = 2 (top lines); the system is therefore in the paramagnetic phase since the 
Griffiths is found for T < T c (l) ?s 1.8205. Time is measured in updates per spin. 
Solid lines represent results of the RS theory. Dashed and dotted lines denote 
the averages and averages ± standard deviation, respectively, as measured over 
20 MC simulations of systems with TV = 10 6 spins. For clarity we plot only the 
average MC magnetization. The size of the symbols is smaller than the error bars. 
Right: histograms (RS theory) of the two field distributions P(±l, n) measured at 
t = 20 compared to the corresponding MC results (markers with error bars). The 
top and bottom panels refer to the temperatures T = 3 and T = 2, respectively. 



depicted and compared in figures [JO In addition we also compare in these figures 
the theoretical predictions for the histograms of fields Pt(s, n) with the corresponding 
MC results, as measured in the final stage of each simulation. We observe that the 
theory correctly predicts both the trajectories of the macroscopic observables and 
the distributions of fields obtained in the MC simulations. Furthermore, one clearly 
notices the profound differences between the macroscopic dynamics of the model (|49|) 
in the paramagnetic phases (figure [T]) versus the Griffiths phase (figure [3J . 

The mesoscopic picture usually put forward to understand the dynamics of spin 
systems in the Griffiths phase [32] is that of local spin clusters that can be regarded 
as independent from (or only weakly dependent on) the rest of the system. Each 
such cluster behaves as a finite (size n) local ferromagnet, with its own 'local' ordering 
temperature T n . A cluster of size n is more likely to be found in the disordered m„ = 
state (where m n is its magnetization) above T„, and in an ordered m n ^ state for 
T < T n . At low temperatures the cluster is equally likely to be in one of its two ground 
states ±m„, which are related by the reversal o~i — » — o~i of all spins in the cluster. 
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t n n 

Figure 2. Left: evolution of the magnetization and energy per spin for k = 3, 
p = 0.2, J = 1 and 8 = 0. The temperatures are T = 1.5 (bottom lines) and 
T = 1 (top lines), so the system has entered the Griffiths phase. Time is measured 
in updates per spin. Solid lines represent results of the RS theory. Dashed and 
dotted lines denote the averages and averages ± standard deviation, respectively, 
as measured over 20 MC simulations of systems with N = 10 6 spins. For clarity we 
plot only the average MC magnetization. The size of symbols is smaller than the 
error bars. Right: histograms (RS theory) of the two field distributions P(±l,n) 
measured at t = 50 compared to the corresponding MC results (markers with 
error bars). The top and bottom panels refer to the temperatures T = 1.5 and 
T = 1 respectively. 

In order to go from m n to — m n the cluster has to overcome an energy barrier E n . 
The microscopic time r„ required for this operation to occur is given by the Arrhcnius 
form r n ~ exp[—E n /T]. The collective behavior of these clusters is thought to be 
responsible for the slowing down of the dynamics in the Griffiths phase 32J. 

The above picture indeed allows us to interpret the results of the present study. 
Our numerical results (figures [l][3]) refer to regular random graphs of degree k = 3, 
with dilution strength p = i, which is below the percolation threshold p c = \ for 
this graph. The simulated system therefore consists of independent clusters of finite 
size, and the density W n (p) of large clusters decays exponentially according to ((501) . 
The Griffiths phase of the model (for k — 3, J = 1 and p = -|) is the range of 
temperatures < T < T c (l), where T c (l) = 1.820478(6) is the critical temperature of 
the corresponding 'clean' undiluted system. Above T c (l) all clusters are paramagnetic, 
and the magnetization and the energy both relax quickly to their equilibrium values 
m = and (|53|) . respectively (see figureQ]). The distribution of fields P(s, n) (see figure 
[!} is symmetric, i.e. P(s,n) = P(—s, —n), as it should be in equilibrium when 6 = 
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Figure 3. Left: evolution of the magnetization and energy per spin for k = 3, 
p = 0.2, J = 1 and 8 = 0. The temperatures are T = 0.5 (bottom lines) and 
T = 0.25 (top lines), so we have entered further into the Griffiths phase. Time 
is measured in updates per spin. Solid lines represent results of the RS theory. 
Dashed and dotted lines denote the averages and averages ± standard deviation, 
respectively, as measured over 20 MC simulations of systems with N = 10 6 spins. 
Right: histograms (RS theory) of the two field distributions P(±l,n), measured 
at t = 100, compared to the corresponding MC results (markers with error bars). 
The top and bottom panels refer to the temperatures T = 0.5 and T = 0.25 
respectively. 



in In the Griffiths phase, in contrast, both paramagnetic and ferromagnetic 

clusters are present. For short times the paramagnetic and ferromagnetic clusters 
evolve to the m n = and m n ^ states, respectively. At intermediate times the 
magnetizations of paramagnetic clusters simply fluctuate around m n = 0, whereas the 
ferromagnetic clusters will 'flip' m„ — > —m n as frequently as the relaxation time of 
the cluster r n allows. Larger clusters require more time to 'flip' due to their energy 
barriers being proportional to their sizes. Furthermore, for lower temperatures the 
ferromagnetic clusters tend to stay longer in each of their two ground states ±m ra . 
Eventually the whole system ends up in the zero global magnetization state. However, 
how quickly this would happen depends on the control parameters of the system. In 
the high temperature region of the Griffiths phase the ferromagnetic clusters would 
'flip' frequently, and although the relaxation time of the order parameters decreases 
at lower temperatures, it is still relatively quick; see figure [U We also observe in this 
latter figure that the energy attains its equilibrium value, given by (|53p , much earlier 
than the magnetization. This marks the onset of the main stage of the dynamics, 
where only the flips m n <-> — m n of the ferromagnetic clusters are relevant. As we 
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Figure 4. Evolution in time of the energy per spin E and the magnetization 
m, now shown as trajectories in the (m, E) plane, for k = 3, p = 0.2, J = 1, 
8 = and temperatures T = 1,0.5,0.4,0.25 (from bottom to top), all of which 
correspond to the Griffiths phase. Solid lines represent the predictions of the RS 
theory. Dashed lines denote average values measured over 20 MC simulations of 
systems with N = 10 8 spins each. The simulations were run for WON sequential 
spin updates, for all temperatures, and the theoretical predictions calculated for 
the equivalent real-time duration t £ [0, 100]. 

decrease the temperature further the dynamics becomes very slow, see figure [3] Here 
the energy has attained its equilibrium value, but the magnetization has not. The 
number of ferromagnetic clusters has increased, and so have the relaxation times t„ of 
those clusters which were already ferromagnetic at higher temperatures. Furthermore, 
at T — 0.25 we observe that in the MC simulation the equilibration times diverge 
with the system size TV (at T = 0.5, in contrast, the system can be still equilibrated 
on timescales significantly less than the system size). This suggests the existence 
of another critical temperature T*, which for the parameters of the system in this 
study would be located somewhere in the interval 0.5 < T* < 0.25, that separates 
the Griffiths phase into two further distinct regions of relatively slow and relatively 
fast dynamics, respectively. A possible mechanism behind this further (dynamic) 
transition would be that the number of clusters which are ferromagnetic becomes 
extensive, combined with diverging cluster relaxation times. Interestingly, the flow in 
the energy versus magnetization plane (see figure [J) for the temperature T = 0.25 
is distinct from that observed at higher temperatures, in terms of an apparently 
discontinuous direction change, and this is observed in both theory and simulation. 
In contrast, in the paramagnetic and high temperature Griffiths regions of this model 
the trajectories in the (m, E) plane are smooth. It is not yet clear to what extent 
the temperature at which the distinct direction of (m, E) flow sets in is related to the 
suggested dynamic transition temperature T* of diverging relaxation times. 
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5. Summary and conclusions 

In this paper we built on a recent study |17j in which a dynamical replica theory (DRT) 
was developed to solve the (sequential) stochastic dynamics of finitely connected 
Ising spin systems with random bonds. Here we generalized this theory to include 
systems on random graphs defined by arbitrary vertex degree distributions (as apposed 
to the Poissonnian ones of [17] ). We have used the exact dynamical equation for 
the joint spin-field probability distribution, that was derived in [17] ■ as a starting 
point. We closed this equation following the standard assumptions of DRT. The 
resulting macroscopic theory takes the form of a nonlinear diffusion equation coupled 
to a functional saddle-point problem, where the latter involves replica density order 
parameters that are to be solved at each instance of time. We showed that the results of 
equilibrium statistical mechanics [34] can be recovered within our dynamic theory, and 
that the equilibrium solution of the model is a stationary point of our macroscopic 
equations. The saddle-point equations resulting from making a replica-symmetric 
ansatz can be solved numerically by a population dynamics algorithm [7] . The results 
in [17j for random graphs with Poissonian degree distributions are easily recovered 
from our generalized equations. 

We have applied our theory to the dynamics of the diluted Ising ferromagnet in 
the Griffiths phase. This model is an Ising ferromagnet defined on a random regular 
graph from which edges are removed randomly, with some probability 1 — p. The 
local fields in this model take integer values, which simplified our dynamic theory 
to a system of ordinary differential equations for the joint probability distribution 
of Ising spins and integer fields. The functional order parameter of the saddle-point 
problem is a distribution over real- valued 2x2 matrices. We have solved our dynamic 
equations numerically for random regular graphs of degree k — 3 with dilution p = ^ , 
and calculated the evolution in time of the magnetization and the energy per spin 
in both the paramagnetic and the Griffiths phases of this model. Dynamic Griffiths 
effects are clearly present in the Griffiths phase. The magnetization equilibrates much 
slower than the energy, and this discrepancy becomes even more severe in the low 
temperature region of the Griffiths phase. In contrast to the paramagnetic phase and 
higher temperature region of the Griffiths phase, the energy per spin appears to be no 
longer a smooth function of the magnetization in the low temperature region of the 
Griffiths phase. The equilibration times of the MC simulation, the results of which are 
in good agreement with the numeric solutions of our theory, diverge with the system 
size in the low temperature region of the Griffiths phase. 

The predictions of the dynamical theory presented in this paper for the diluted 
Ising ferromagnet in its Griffiths phase are remarkably accurate. To us this is not 
entirely surprising, for at least two reasons. First, the dynamic replica theory and its 
variations have in the past already proven to be very accurate for ferromagnetic models 
on regular [12l [14] and Poissonian [17] random graphs. Second, the extended version 
of DRT considered here describes the evolution of the joint spin-field probability 
distribution. In equilibrium, the cluster expansion of the magnetization derived in 
[22] can be recovered within the cavity approach [37] . which is equivalent to the 
replica method. This suggests that the joint spin-field probability indeed contains 
the relevant information about the clusters which is responsible for the slow dynamics 
in the Griffiths phase. 
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Appendix A. Averaging over disorder 



In this section we compute averages over disorder {cijJij} in the equation (p~5|) . 
First, we exploit the i <-> j symmetry of the interactions CijJij to write the disorder 
dependent term of this equation in form more convenient for further manipulations. 
Second, we write the Kronecker delta in the definition of the connectivity {cy } disorder 
(O in its integral representation (fT5|) . This gives us 



(• ■ -){c I3 J l3 } 



cu'5{h — H\i + 2 Jui s] e 



{CijJij} 



C i<j 



N' 



(A.l) 



Taking the average over connectivity disorder {Cy} leads us to 



n 



dk 



1 -iey . h 
— e ^ a >* 



c«/ 5 [h — H}, + 2 J w s] e 
dfc 



En +(i-^,o 

J C z<j 

>E,<j c a [■JijY.c 



{■/«} 



n 



2?r 



x — / dJ P(J) - ffj> + 2Js] e 



><n 



dJ P(J) e~ iJ ^« {*?"?+&?"?}-Hfe+**} + (l--^)j (A.2) 



where in the last line of above expression i ^ I and j ^ Finally, upon re- 
cxponcntiating (|A.2|) we obtain our desired result for the disorder average in (p~3|) : 



(■ ■ •){*jJij} 



n 



dfc 

2^ 



1 gife.fci 



x — / dJ P(J) <5[/i - + 2Js] e 



x exp 



c 
2~/V 



{h?o-°,+?i°,(Tf}-i{fc«+fc € /} 

^ |y*d jp( j) e - ij ^° } -ufci+fcj ■> _ xj + O (jvo) 



(A.3) 



Dynamics in the Griffiths phase of the diluted Ising ferromagnet 
Appendix B. Computation of the kernel A[. . .} 



21 



In this appendix we give details of the calculation which leads to the path integral 
(fT9jl . We insert our result for the term (cwA u >[. . •])/ Ci .j ( .i (O into the sum (fTU|) , 
which gives 



A[s, s ;h,h ;s 



Z N 2 



x exp 



E/n 

i^£> J rha 



d£>q(T, h) 
2ir/N 
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\cN [ dhdh'dkdk'P(cr, h, k)P{<r', h, k')x... 

-U[h.(T'+h' .a]-i[k+k'] - 1^ + 0(N°) 

E-E/rT{^}£n[f •** 



x exp [jy^fej. {Hi - 6} 

i 
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Next we rescale the conjugate integration variables according to D a — > D a Ah, and 
define the function 

Ztt 



x e 



- i T, T .h, a AkB«(T,/i) S T ^S[h-H"]+ik i ki-iP{(Ti,h i ,k i ) 



(B.2) 



insertion of which into (|B.1[) . followed by further manipulations, leads us to 

1 r 1 i N r - - 
A[s, s'; h, ti; s] = - [—J / {dPdPdD} 



x exp 



iN Y AhD «(T, h)D{r, h)+iN^J dhdkP{cr 7 h, k)P(a, h, k) 

-cN Y [ dhdhdkdk'P(cr, h, k)P(tr', h, k')x ... 
2 cr,(T> J 

x ^ c -u[h.cr'+h' .cr]-i[k+k'] _ ^ 
E lQ gE / dH ^ dh ^ f &i M[H l ,h l ,a l \k t ,k t ,9} + 0(N°) 

IT J J— IT 
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,<Te\ki,k e ,9] 

t> M [J2> ,h v ,a v \ hf , %i> , 9} 

-ij[ht.tr t i+hir .<Tt] 



Y [dHz.dhz, f dh> 



x 5 s ,^ l 5 s ^x5[h! - H}] (6[h-H}, + 2J~s] e -vin*-<rt> 

[kt+k e ,] 

dHfdhf 



x e 

x 



x Yl I dH e ,dh 



dkgM[He, he, cre\ke, kt, i 



dki> M [H e > ,hi>,trt'\ kt> , ki> , 9} 



(B.3) 



Now the terms in the sums over i, £' variables are dependent only on the random 
connectivity variables {ke,k'g\, which are independent and distributed according to 
P c (k), hence by the law of large numbers we arrive at the result (fl9|) . 

Appendix C. Calculation of the Fourier transforms 

Here we compute the Fourier transforms J dh e~ lX -h dk e~ lkm M[H, h, cr\k, k 



of the function M defined in equation ([ZTjl . where x £ K ra and m G Z. First we 
expand that part of the exponential function which depends on Q, which is defined in 
equation (|26|) . Next we integrate out the k variables, which leads us to 



dh e" 



iX.fl 



dk e- ikm M\H,K(r\k,k. 



= [ dh jh.{H-x-e}-iY, a b a {a a ,H a 



Q k ' m (cr,h) (C.l) 

(k — my. 

In the above we used the short-hand Q(cr, h) = Q(cr, h, 0) + 1. Raising Q to the power 
k — m gives 

Q k - m (<r,h) 



<7' 

k—m 



Y / dhdk'P(a',h,k') e~ ik ' ( e" 1 



\j\h.cr'+h .er] 



n 



Y / dh e dJ e P(J e 
ih. Y^e J i0 1 



&k t P(at,ht,kt) e~ ik * e- u * h * a 



(C.2) 

Now inserting above result into the expression (|C.1|) and integrating out the h variables 
yields equation for the Fourier transform (l30|) . 

Appendix D. The joint spin-field probability distributions in equilibrium 



In this section we compute the joint spin-field probability distributions D and A in 
equilibrium. We note that both can defined via the Fourier transforms (f30| of the 
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function M. First, we consider equation for D (f2"51) . Using expression (|28p for the 
conjugate parameter D a in this equation, combined with the equality (|33p . gives us 

k 

D(a, /*) = E p d k )^ E f dH LI E / dhtdJ e P(J e ) 



fc>0 



x 8 



l L er t 



dk e P{<Te,h e ,k)e- i ~ k 'e- u * fl ' (T 



E^E/^n 



fe>0 



t=i Lcr 



J2 / dJ e P(J e 



dk e P(<rt,k) e~ ik 

r 

(D.l) 



where M k is defined in (|32j) . Summing and integrating over the variables a 1 and H 1 , 
respectively, leads us to the equilibrium form ([36]) of the joint spin-field distribution. 
In a similar manner we obtain the equilibrium version of A, which is given by 

P c (k)k P c (k')k' 1 



A[s,s';h,h';s] = fdJP(J)^-^ 

J A k,k' 



M k M k 



' rr rri J 



k-1 



n 

fe'-i 
n 



E dJ t P{Ji) / dktPfakde-** 



dJ r P{J r ) / dk r P(er r ,h r )e-' lK 

x e 0(T.H+f3<T'.H'-pjCT.CT' 

x 5s', a t S s>< 8[ti - Hi] 5[A - ffi + 2Js]. 
The above result can be written in the form given by equation (f37l 



H-^Jucri-O-Ja' 

i 

H l — ^ y J r o ' T — — Jcr 



(D.2) 



Appendix E. Stationary points of the dynamic equation 

Here we show that the probability distributions D and A in equilibrium are stationary 
points of our dynamic equation |8]). First, we consider that part of ((8]) which is 
dependent on the joint spin-field distribution D(s, h) only. Inserting the equilibrium 
form (f56"| of this distribution into the first line in the right-hand side of ([5]) leads to 

i[l + stanh[/3/i]] e- 0ah ^{h]-hl-sta,nh[/3h}} e 0sh ^[h] 

= ®[h]{- sinh [f3s h] + cosh[f3sh] tanh[f3sh]} 

= 0. (E.l) 

Second, we compute that part of the right-hand side of (jHJ) which is explicitly 
dependent on the kernel A[s, s', h, h', s) only. Using our equilibrium form ([57)1 of this 
kernel in the last two lines of the right-hand side of ([8]) results in 

[ dti[l - s'tanh[/3/i']]<J e f}sh+l3s ' h ' +fiJss ' K[h + Js'-h' - J a]} 
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~^cJ2j dh 'l l ~ s ' timh[(3ti}]( e^h+ps'h'-pjss' k][h _ Jg ,. h > _ Jfj ^ 

s' 

= ice' 35 ' 1 JdJP(J) fdh'{ 

+ ([l-tanh[f3h']] C 3 ' 1 ' - [l + tanh[/3/i']] e^') e 0Js A[h + J;h'~ Js] 

'[l+tanh[/3h']] e- ph ' - [l-tanh[/?/i']] e ph '\ e~^ j3 A[h - J; h 
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'-Js]} 



= 0. (E.2) 

We conclude that the right-hand side of the dynamic equation ([5]) is exactly zero for 
all s € {— 1, 1} and all h £ R as soon as the equilibrium relations ([36]) and (|37|) hold. 

Appendix F. RS calculations 

In this section we derive an equation for the functional distribution (|4ip and compute 
the replica symmetric versions of the kernels A and D. First, we compute the 
functional distribution VF[{P}; \m], where m £ Z. For this we consider the Fourier 
transform J^dke~ lkm Pjts{ (T , h, k) of the RS order parameter function. Using result 
(|40[) for Mas, and the saddle-point equation we have 



{dP} / dk W[{P}; k] e~ ikm J] P(a a , h a ) 

J-* a=l 

J dH j% dk e~ ikm M RS [H, h, a\k, k, 6} 



k>0 



Ea/™ M RS [H,h,cr\k,\ 

fe-T 



fc>0 



(k— m)\ 



,J_ 

Ml 



n 



dJ^P(J £ ){dP £ } / dk e W[{Pt};kt]e 



i /»• , 



n „ k — rn „ 

<*=!"' 1=1 erf J 



Jt[hfU a +h a af] 



{dP}^P c (fc) 



fc! 



k>Q 



(k-m)\ 



1 erf 
k—m r 



n 



d.J e P(J e ) {dP,} / dk e W[{P e }; k t ]e 



- i k t 



a. 1 1 



P(a,h)- 



Z[{Px,...,P k ^ m }] 



x 



1 " 
—Z[{P 1 ,...,P k ^ m }] n l[P(a a: h a ) 



(F.l) 



a—l 

where we have used the short-hands 

k r 



n 



dJiPiJi) {dP e } / dfc £ PF[{P £ }; A/] e 



J2 J dHdh ^,ff)e i/l(H - 9) n(E J d hiPfXvi,he)e- Uilile 



(T-\-h(Ti\ 



(F.2) 
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Z[{Pi, . . . , P k -m}} = E / dHd hd(°, H)e ih ( H -^ 
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k—m 



£=1 <Ti ^ 

k — m r 

^EII [E/ d ^^'^) e ~ 1 * / ^ ff ] d ( <7 'E J ^ + ^)- ( F - 3 ) 



<y 1=1 at 



Solving equation (jF.lj) for the functional distribution J* n dk W[{P}; k] e~ lkm , followed 
by taking the replica limit n — > in the functions and Z™ of the resulting 
expression, then leads to equation (|42f . 

Second, we compute the RS joint spin-field probability distribution D(s, h). Using 
the saddle-point equation (|23|) for this distribution, combined with the result (j40| for 
Afjjs, applied to m = 0, gives us 

D(a h)=TP (k) E<T ^ Mrs[H ' -2 <T|fc ' ° ] 6 ^ 6{h Z H t> 



fe>0 



fe>0 



dJ t P(Jt) {dP t } / dfc<W[{P/}; fei] e - *' 



nLi / dtafo, he)e- iJ * h *° d{<r, h) 5{h -J2 e Jm 



xZ[{P l7 ...,P fc }]" (F.4) 

where the functions and Z[. . .]" are defined by (|F.2[) and (IF.3[) respectively. Taking 
the replica limit in equation (|F.4[) leads us to the result (|4"4"1) . 

Finally, we compute the RS version of the kernel (|22|) . We consider numerator and 
denominator in the average over the vertex connectivities in this equation separately. 
Using equality (140]) for the function Mrs we obtain the numerator 

num = E / dHdH'dhdhM RS [H, h,a\k- 1, 0] M RS [H', h, a'\k' - 1, 9] 

(T ,<J' 

x 5 s> ^S St< 5[h' - H x ] (s[h -H[+ 2JS}e^ h - (T ' +h '-^ 



(fc-1) 



; J 1=1 



dJiP{J e ) {dP e } / dk e W[{P e };k e ] e 



-ikf 



e -c c k'-l > 



(k 



^. n 



dj;P(4) {dQ r } / dk r W[{Q r }; k r ] e-^ 



En 



fc'-l 



E 6h l P i {tTi,ht)e- iJlhia 



n 

r=l 



E / dh r Q r (ar,hr)e- iJ > rC 



d(o-,^2 J e o- e + 9 + J a') 

t 

d{a'^j' r a r + e + J<j) 
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x 5 a >, a 8„y 5[ti - J ^ Jcr '] S i h ~ E J '^ T - - J(T + 2JS ] 



k-1 



En 



l a,cr' 1=1 
k'-l 



/ dh £ P e (a e ,h e )e- iJ 



>hecr 



d(a,^2 j£<j e + + Jo') 



><n 

r=l 

and the denominator 



^2 / dh r Q r (a r ,h r )e- iJ 'r hr ' 7 ' 



n-l\ 



d(a',y~] J' r a r + 6 + Ja) 



>(F.5) 



den= Y I dHdH'dhdhM RS [H,h,a\k,ff\ M RS [H',h ,<r'\k', 



e c c k 

k\ 



n 



En 



dJ e P(J t ) {dP e } / dk e W[{P e };k e ] e~ ik 



^ / dhtP e (c- e ,he)e- iJ * h *° 



e -c c k' r J^ 



k>\ 



k' 



En 



n 



^ / dh r Q r (a r ,h r )e- iJ - h " c 



dJ' r P{J' r ) {dQ r } / dk r W[{Q r };k r ] e" 



(F.6) 



Combining these latter two results in (j2"2"|) and taking the n — » replica limit gives 
equation (|45|) . 



Appendix G. Dynamic equation for the Ising ferromagnet with dilution 

Here we show that the macroscopic equation JSJ for the Ising spin system governed 
by (|49p can be reduced to a system of ordinary differential equations. In the present 
Ising ferromagnet with dilution the fields §5§ can take only discrete values, which 
implies that the distributions iJTJl and (fTOj) can be written in the form (|60|) and (|6T|) 
respectively. Inserting (|60|) and (|6"Tj) into both sides of ((8|) gives 

^ Pt(s,n) 6{h~Jn-8) = ~ [1 + stanh[/3/i]] ^ P t {-s,n) 8(h-Jn-6) 



n= — k 



n=—k 
k 



- [1-stanh^/i]] ^ Pt(*,n) S{h-Jn-6) 



— k 



k-l k-l 



+ \ k E / ^'E 1 - s ' tanh[/3/x']] ]T E 

s' n— — n' — — fe+1 

x (A t [s,s';n,n'|r] Jri-0- Jts] S[h- Jn-6+Jrs']) T 

-lkJ2jdh'[l-s' tanh[/%']] ]T ]T 

s 7 n=— u'= — fc+1 

x (A t [s,s';n,n'\T] 8[ti -Jn' -6-Jts] S[h-Jn-9-JTs']) r (G.l) 
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in which the averages over r refer to the distribution P(t) = pS T i + (1 — p)S T ^. We 
move the time derivative inside the sum on the left of the above equation. On the 
right side we average over r, take the sums over s' , and integrate out h! variables. 
These manipulations produce 

^Pt(s,n) S(h-Jn-6) = ~ [1+s tanh[/3( Jn+0)]] ^ P t (-s,n) 5(h-Jn-9) 

n——k n= — fe 

k 



- [1-s tanh[/3( Jn+0)]] ^ P t (s,n) S(h-Jn-O) 



k-2 fe-1 

+ ^ Y -kp[l-tsaih[l3J(n'+s)+06]]A t [s,lin+l,n'\l]S(h-Jn- 



2 

n— — k n' — — k-\-l 



+ -kp[l+tanh[/3J(n'+s)+06]]A t [s,-l;n-l,ri'\l]6(h-Jn- 

n=-fc+2 n' = -k+l 
k-2 fe-1 ^ 

Y -kp [l+t&nh[f3J(n' + s)+f36]}A t [s,-l;n+l,n'\l}5(h-Jn-6 

n—~k n' = — fe+1 
k fe-1 . 

- ^ ^ -fcp [l-tanh[/3J(n' + s)+/30]]A t [s,l;n-l,n'|l]J(/i-Jn-f 



n=-fc+2 n' = -k+l 

(G.2) 

The result (|64[) follows immediately from the above equation. 
Appendix H. Initial conditions 

In this appendix we compute the relevant initial conditions for the system of ordinary 
equations We choose an initial microscopic state of the system in which each 

spin <Ji is drawn randomly and independently according to Po( cr i) = |(1 + Cimo), 
where mo 6 [—1,1] is the prescribed initial magnetization of the whole system, i.e. 



Po(t)= n^l + o-imo). (H.l) 
j=i 

Given (|H. 1|) . the spin-field probability distribution Pq(s, n) for large Ising ferromagnets 
defined on random graphs with vertex degree distribution P c (k') is given by 

1 N 

P (s,n) = Jim ^PqH-^VA.E^c^ 



i 

i/ 



5 vV n .(H.2) 



= i(l+ S mo)£P c (fc')fl £~(l+a,m ) 

fc'>0 1=1 L o-f 

For the model (|49|) in particular, where is the connectivity of the random regular 
graph and p is the dilution, the vertex degree distribution is binomial 

k 



P 



{k')=( k kl )p k \l-p) k - k '. (H.3) 
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Solving equations (|H.2|) and (fol?)) for the functional distribution W[{P}|1] and the 
function d(s, Jn + 9) then gives 

W[{P}\1] = Y[s[P(a\a')-^(l + am Q ) 



d(s, Jn + 6) = -(1 + smo) 



(H.4) 
(H.5) 



which is the trivial solution of equation ([58 



Appendix I. Population dynamics 

The joint spin-field probability distribution P t (s,n) of the diluted ferromagnet (|4"5)) 
evolves in time according to the system of ordinary differential equations (|64p . Solving 
this system requires computation of the kernel (|63[) . which is dependent on the 
functional distribution W and the function d (the order parameters). The saddle- 
point equations (|57|) and (|62|) establish relations between these parameters and their 
dependence on P t (s,n). However, solving these equations analytically is generally 
ruled out, and one has to solve them numerically using population dynamics [7]. 

The population dynamics algorithm was also used in the preceding version of 
the dynamical replica theory, as developed for Poissonian random graphs [17] . Here, 
however, we take an approach which is slightly different from the one in [17j . We 
note that in our dynamical theory we use P*(s, n) to estimate the order parameters W 
and d. In particular, the values of the order parameters are considered to be 'good' 
when the saddle-point equation (|5T[) for the functional distribution W is satisfied, and 
the probability distribution P(s, n) which is computed via saddle-point equation (|62[) 
equals the instantaneous distribution P«(s, n). This suggests that the change made by 
any numerical algorithm to the order parameters W and d has to reduce the 'distance' 
between the distributions Pt(s, n) and P(s, n), subject to the constraints ([ST]) and (162"]) . 
The Kullback-Leibler (KL) divergence 

p^L(Piip) = EE p *( s ' n ) lo g 

s n 

can play the role of a distance in this context, and we may use e.g. a gradient descent 
algorithm to minimize this distance, viz. 

^^ + e) = - dd{ J n + e) D KLm p) d.2) 

where e defines an 'algorithmic time'. To solve equations (|57|) and (|62j) . we use 
a combination of both population dynamics and gradient descent. To implement 
the population dynamics we create a population of N 2x2 matrices Pi(a\a'), 
where i — 1 . . .A/", and we initialize the function d(s, Jn + 9) for s € { — 1, 1} and 
n G {— k, . . . , k}. The initial values of population {Pi(a\o~')} and function d(s, Jn + 9) 
are set to (|H.4[) and (|H.5[) . respectively, at t = 0. For t > we simply reuse values 
from the previous time step. We then update the population of matrices and the 
numbers d(s, Jn + 9) until they are stationary, via the following process: 

(i) a number k' is drawn from the binomial distribution Bk-i(k') (|59p 

(ii) k' members Pi(a\a') are selected randomly and independently from the 
population 



Pt(s,n) 
P(s,ri) 



(LI) 
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(iii) a new value for P(a\a') is calculated according to 



Pnew(cr\a') = 
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(1.3) 



E." n,=i {E CT1 Pi{*M")} d(a",JT,li *1 

(iv) a member of the population is selected randomly, and replaced with the newly 
computed value P new (a\a') 

(v) a new function d(s,n) is computed according to 

d(s, n) 



,(s, n) — d(s, n) + Ae- 



[P t (s,n) - P(s,n)] 



(1.4) 



T + d(s,n) 2 

where < Ae -C 1, and P(s, n) is computed according to ([62]) by averaging over 
the instantaneous values of the population. 

The rule (|I.4p used to update d(s, n) can be regarded as an approximation of the 
gradient descent equation (|I.2|1 . which can be derived as follows. First we use the 
definition of the KL divergence (II. 1[) and equation (f6"2")l for P(s, n) to compute the 
partial derivative in (|I.2|) . giving (with the short-hand d(s, Jn + 0) — > d(s,n)) 

9 D KL (P t \\P) = ~ PtM 



dd(s, n) 



d(s, n) 

o'to' ^ ' Ti Tl. 



,r fe ) /n[{dA}w[{A}|l] 



(=1 a i 



K 



(1.5) 



The result p.5[) takes a very simple form when there is no disorder and the distribution 
W[{P}|1] is a functional delta, where one is led to 



-r-d{s,n) = — 

de a{s, n) 



[P t (s,n)-P(s,n)]. 



(1.6) 



To reduce computational costs we use in our population dynamics algorithm 
approximation (|I.6|1 . rather than the full version of the gradient descent p.2[) which 
would have required computation of (|I.5p . First, however, expression p.6[) is slightly 
modified according to l/d(s,n) — > d(s,n)/[l + d(s,n) 2 ], to prevent unbounded 
increasing (or decreasing) of Ae in the discrete version of (|I.6[) . The number of 
iterations required to solve saddle-point equations (|57I62I) by the algorithm presented 
in this section was found to be typically of order 10 2 7V, for the population size M = 10 4 . 
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